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Abstract 

A new method is presented for the computation of electromagnetic scattering from 
axially symmetric bodies. To allow the simulation of inhomogeneous cross sections, the 
method combines the finite element and boundary element techniques. Interior to a 
fictitious surface enclosing the scattering body, the finite element method is used which 
results in a sparse submatrix, whereas along the enclosure the Stratton- Chu integral 
equation is enforced. By choosing the fictitious enclosure to be a right circular cylinder, 
most of the resulting boundary integrals are convolutional and may therefore be evaluated 
via the FFT with the system is iteratively solved. In view of the sparse matrix associated 
with the interior fields, this reduces the storage requirement of the entire system to 0(N ) 
making the method attractive for large scale computations. This report describes the 
details of the corresponding formulation and its numerical implementation. 
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Chapter 1 

Introduction 


A restraining factor in the numerical simulation of three-dimensional structures for 
electromagnetic scattering computations is the storage requirement associated with the 
chosen method. For sub- wavelength structures traditional methods [1] have been found 
to work well. However, for structures spanning several wavelengths, the storage require- 
ment limits the use of these methods. 

For the special case of axially symmetric structures or bodies of revolution (BOR), a 
reduction of the storage requirement is accomplished by reducing the three-dimensional 
problem to a set of two-dimensional ones. Several moment method codes have been 
developed for the solution of these ([2] - [7] and others). However, for large structures 
the required storage of 0(JV 2 ), where N denotes the number of unknowns over the BOR 
cross section, limits their use. 

To further reduce the storage requirement, hybrid finite element methods ([8]-[12], 
etc.) may be used, since the storage associated with the finite element method is 0(N) 
in contrast to the 0{N 2 ) requirement of moment methods. These methods differ from 
one another primarily by the application of the radiation condition. The most accurate 
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method enforces the radiation implicitly through an application of the boundary integral 
equation over a fictitious boundary enclosure [11], and in this case the storage is still 
0(JVj ), where Nt is the number of unknowns on the boundary. However, through a 
judicious choice of the enclosing boundary, the storage requirement can be reduced to 
O(N). This can be achieved by selecting the enclosing boundary to be rectangular or 
circular [15], [16], making some of the integrals convolutions which can then be evaluated 
via the FFT when an iterative solution scheme is employed. 

The proposed method combines the Unite element and boundary element methods for 
the solution of inhomogeneous bodies of revolution. The coupled potential equations [10] 
are discretized via the usual finite element method, and the resulting system is augmented 
by a discrete representation (via the boundary element method [13]) of the Stratton-Chu 
equations [14]. By choosing a right circular cylinder to enclose the scatterer, some of 
the integrals become convolutions and their discrete counterparts are then evaluated 
via the FFT in conjunction with an iterative solution procedure as was done in the 
two-dimensional case [15]. With some care, the storage is reduced to O(N). 

In this report, we describe the formulation for the proposed finite element - boundary 
element method as applied to the body of revolution. Some preliminary results are shown 
to be in reasonable agreement with the method of moments (MOM). 
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Chapter 2 

Analysis 


Consider the body of revolution (BOR) illustrated in fig. 2.1. To employ the proposed 
finite element - boundary element (FE/BE) method, the BOR i6 tightly enclosed in a 
fictitious finite length cylinder, which divides the entire space into two regions, i.e. the 
one enclosed by the cylinder and the other exterior to it. Since the interior region is 
generally inhomogeneous, the finite element method is suited for formulating the fields 
of that region, whereas the boundary element method is applicable for the exterior free 
space region. A usual approach [3] for treating BORs is to introduce a Fourier series (in 
the azimuthal coordinate <j > ) representation of the fields, reducing the problem to a set 
of two* dimensional ones. The finite element - boundary element method is then used to 
compute each modal field and the final result is found by adding the modal fields. 

In the following, we present the finite element and boundary element formulations 
for each mode. First, the finite element formulation is developed. 
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Figure 2.1: General surface of revolution. 
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2.1 Finite Element Formulation 


In this section, we derive the analytical coupled azimuth potential (CAP) equations 
[17] which are then discretized via the finite element method. 

2.1.1 Analytic CAP Formulation 

Maxwell’s equations in a source free region (a time dependence is assumed and 
suppressed) are given by 


V x E{r) — —jufiH 

(2.1) 

V x ~B(J) = jueE 

(2.2) 

V • ~D(r) = 0 

(2.3) 

v-%(7) = 0 

(2.4) 


(2.5) 

For axially symmetric media, the fields may be represented as 

Fourier series in the 

cylindrical coordinate <f> as 


E(?)= £ mp,*)^"* 

m*— oo 

(2.6) 

n%( 7 ) = 53 ^(p,*)^” 1 * 

m=— oo 

(2.7) 

and when these are substituted into Maxwell’s equations (2.1) and (2.2), we obtain 

£ \jmh mz - &(Rhmt)] = j€ r e mp 

(2.8) 


(2.9) 


(2.10) 
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jj = “j f r c mz 


( 2 . 11 ) 

( 2 . 12 ) 


i [j Wm^ — ^(-RCm^)] — JMr^ro* 


(2.13) 


with 


R = kop (2.14) 

Z = *o* (2.15) 

to be referred to as normalized coordinates. Substituting hm X of (2.13) into (2.8) gives 

c mp — jfm [TO ^j(-^ c m^) 4" Rf*r ^'(•^bn^)] (2*16) 

where 

fm = [.RV-m 2 ]" 1 (2.17) 

« 2 = /z r c r (2.18) 

Substituting e mp of (2.8) into (2.13) we obtain 

frm* = jfm [rn &(Rh m4> ) + Re r &(Re m4 ,j\ (2.19) 

Substituting h mfi of (2.9) into (2.12) yields 

Cm* = jfm [to &(Re m4) ) - Rp r ^(Rh^)] (2.20) 

Substituting e mj of (2.12) into (2.9) yields 

h-mp — jfm [to •£p(Rh m 4 > ) — Re r •£g{Re m 4 >j] (2-21) 
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Equations (2.16) through (2.21) may be written in compact form as 


x ?(R, Z ) = jf m \m4> X V«V»e - HrKVith] (2.22) 

4 > X H(R, Z) = j/ m x + CrfiViV»e] (2.23) 

^•*0R,Z) = ^ e /-R (2.24) 

4>-H(R,Z) = fa/R (2.25) 

where 


^ A + * & (2.26) 

Rewriting (2.10) and (2.11) as 

■EV* • {4> x 7^) = -y'€ r ^e (2.27) 

RV t • (^ x e m ) = j> r V^ (2.28) 

and then substituting (2.22) and (2.23) into them, we obtain the CAP equations 

V t • [/m(fr-RV t V> e + mj> x V t ^ fc ] + = 0 (2.29) 

Vt • [/m(Mr-RV t V»A - X V t Ve] + “p = 0. (2.30) 

This system may be written in operator form as 


where 


and 


Lj> = 0 


v t . [f m c r RVt] +£ mV,[/ m £ X V t ] 

x V,] V, • [f m fi r RV t ] + $ 


(2.31) 


(2.32) 


V> = [V>e V’/i] T 


(2.33) 
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Figure 2.2: Cross section of a generating surface enclosed by the fictitious boundary r„. 
2.1.2 Discretization of the CAP Equations 

To discretize (2.31), we first enclose the generating contours of BOR in a fictitious 
boundary r o and the axis of symmetry as shown in Fig. 2.2. The contour T 0 is chosen to 
be rectangular in shape thu6 generating a right circular cylinder. The region interior to r o 
is divided into N e linear triangular elements and within each element the corresponding 
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weighted residual expression is written 


// 


RN?(R,Z)Rn'dS e = 0 


(2.34) 


where ( RN f) is the weighting function and Rfa is the residual. Further, Nf is the usual 
linear shape function found in any finite element book [18]. Using this definition, (2.29) 
and (2.30) may be written 



and upon using the identity 


(2.35) 

(2.36) 


(RN')V t ■ A 3 = V, • (RN*A 3 ) - A* • V t (RN?) 


(2.37) 


we obtain 


(2.38) 


(2.39) 


/ f RN f [V t • {RN?f m URV^ + Wx V t i > h ) } 

S e 

+(r1> e N? - f m (e r RV t rp e + mj>x V^a) • V t (J2JV?)] dS* = 0 

JJ RN f [v, • {RN?f m {^RV t rj> h — m4> x V,*.)} 

Sr 

+firil> h N? - f m (fi r RV t tl> k — m$ x V,V> e ) * V t (RN?j\ dS' = 0 
Further, by invoking the divergence theorem (2.38) and (2.39) may be written as 

JJ { [-fm (cr RVtrpe + x V,^)] . VtiRN?) + €r^Nt} dS e 

S e 

+ J C '*‘ [RN'f m (trRVt^t + m<£ x V«V>fc)] dl e = 0 (2.40) 

//{[-'» (fi r RV t if>h - m4> x V«ifr e )] • V t (RN?) + /i r ^ fc AT?| dS e 

+ J ft • [-R-W,' / m (n r RVtll>h - m4> X VtV'e)! d/ e = 0 (2.41) 


5 * 
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where ft is the outward normal along the boundary, C' y of the eth dement. Finally, 
these may be simplified by making use of (2.22) and (2.23), yidding 


// { f-/ m + ™j>X V,**)] * V,(JW?) + t^Nf) dS' 

S' 

-l RNftthmtW = 0 (2.42) 

Jc • 

//{[-A. (flrRVt'I’h -MX V,*.)] • V,(^) + MrtW} <*£' 

5* 

+ i RN?{je mt )dl' = 0 (2.43) 

./c* 

where 

e m < = i ■ e m (2-44) 

h m t = i-h m (2-45) 

with 

t = hx (2-46) 


To form a system of equations over the eth dement, the fidds are represented as a 
linear basis expansion as 


iJ> e (R,Z) = '£4>: j RNf(R,Z) 

3 - 1 

MR,z) = '£'i’h RN !( R ’ z ) 

;=i 

Substituting these into (2.40) and (2.41) yidds 


E 


3=1 


I! { [-fm<rRV t (RNn . V t (RN‘) + <rN‘N'\ 1>' ej 
S' 

—Tnj> x Vt(RNf) • Vt(RNnP hj } 


(2.47) 

(2.48) 
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(2.49) 


3 

E 

Ly* 


~<f RN?(jhmtW = 0 

Jc* 

JJ { [-f m HrRV t (RNn • V«(JW/) + HrN?N'\ 

+m4> x Vt(RNj) • V t (RN t ni>lj} dS*] 

+ l RN?(je mt )dl c = 0 
./c* 


(2.50) 


lb proceed further, it is necessary that we evaluate the integrals over the surface of the 
element. 

Assuming e r and fi r are constant within a given element, (2.49) and (2.50) 

E [<«&<( - - £. XSfdhmiW = o (2.51) 

E[W + «/jW RN’Uc^JF = 0 (2.52) 

jml JC ' 

where 


*?, * JJ \-fmRV t {RNf) . Vt(RNj) + N*Nj] dS e 

S e 

% = JJ [W X Vt(RNf) ■ V t (RNnr hj ] dS' 


(2.53) 

(2.54) 


Summing over all elements to obtain a solution for the entire problem region, our system 
becomes 

N e 3 „ N a 


1C S [ £ r a 0'^e; “ — 0 

c=l ,7=1 *=1 

5Z 5Z [ 6 0^ei + ria'j'Ph,] + 51 = 0 


(2.55) 

(2.56) 


e=lj=l 


#=1 


where 



and i? is the pulse function equal to unity in the sth element. Note that in (2.55) and 
(2.56) the contour integral contribution canceled out except along the boundary T 0 as 
shown in Appendix D. 


In block matrix form (2.55) and (2.56) may be written 


Ka 

*al 

Ah 

Ah 

Baa 

-Bal 

-Bax 

-Bad 

0 

—Caa 




Ah 

Ah 

A e id 

-Bu 

-Bu 

-Bu 

-Bid 

0 

0 



A',a 

Ah 

Ah 

Ah 

-B xa 

-Bxi 

-Bxx 

-Bxd 

0 

0 



Ah 

Adi 

Ah 

Ah 

Bda 

-Bdi 

-Bdx 

-Bdd 

0 

0 



Baa 

Bal 

Bar 

Bad 

Ah 

Ah 

Ah 

Ah 

Caa 

0 


kmfa 

Bla 

Bu 

Bu 

Bid 

Ah 

Ah 

Ah 

Ah 

0 

0 



B, a 

B.i 

B xz 

Bxd 

Ah 

Ah 

Ah 

Ah 

0 

0 


^m^z 

Bda 

Bdi 

Bd» 

Bdd 

Ah 

Ah 

Ah 

Ah 

0 

0 














jCmt 










. 


jhmt 







-fo 

0 0 

0 0 

0 0 

0 

0 0 


] (2.58) 


in which we have substituted V'e and V’fc with e m ^ and h m ^, respectively, and 


N, 



(2.59) 

eel 


N. 



(2.60) 

eel 


N t 


b ='£*•;, 

(2.61) 

e=l 


N* 


'<? 

Wi- 

ll 

o 

(2.62) 
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The subscripts on A t,M , B and C refer to the various regions of JJ and its boundary. The 
elements of (2.59) • (2.62) are derived in Appendix E and are listed as follows 


«'i = - {?!<*) + Fj**1)Qn - 2(7? a i + 7 J ? «?)Qao - 2 + fitf)Q 2 1 

-PtFjQn - (4 7 J h] + ff/?)$30 

+Q?a^Pio + (tfa* + ^o?)P u + ( 7 * a* + 7* a ? )^ao + (#7; + ^7i)P 2 i + A'/JJPia 

+7? 7i P 3°]( 2 ft e )2^ 2 - 63 ) 

and 


m 


IW -««;»»+ )«*>] (2.64) 

where the Ps and Qs are defined in Appendix E- The elements of C are 


C„ = cJi 
£«+!,* = C$1 


(2.65) 


where 


«5i = * = ^ 

Ci, -i(A- ± R[)(R\ + Ri) =fs(K? + r;r; + R?) 

= *[**;(«; + as) -h*; j + *; Ri + «j J )i 

To form a complete system, (2.58) mu6t be appended by a discrete version of the 
boundary integral equation to be discussed next. 


for r a2 

upper sign for T 0 i (2.66) 
lower sign for 
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2.2 Boundary Integral Formulation 


The electric and magnetic fields are represented in the unbounded region by 

T(r) = T'(t) + T(t) (2.67) 

~S(7) = H\r) + TP (7) (2.68) 

where £* (r) and TT (r) are the incident fields and the scattered fields are given by the 
Stratton-Chu equations [14] 

= ff {-jwn [«' X F(r / )] g{7,7 > ) + [»' • E(r / )] V'g(r,r') 

S' 

+ |n' x T(F0] x dS' (2.69) 

Z'(7) = ff {M [«' X TCr 7 )] $(?,?') + [n' • ¥(f / )] V'g{Tf) 

S' 

+ [n' x ¥(r / )] x V'ff^.f')} dS' (2.70) 


where 7* and 7 are the source and observation points, respectively and 

-iJkoF-T’l 

^ (2J1) 

is the free space Green’s function. It is convenient for computational purposes to elimi- 
nate the presence of the normal held components and after 6ome manipulation we obtain 

= ff {—>*» [*' X g(T,r) + J- [«' . V' X vy?,?*) 

s 

+ [n' x E^)] x V' 5 (r,r')}dS' (2.72) 
**>E» = ff fa [»' x Tff)] g(7,Y) -i[n'-Vx EfT)] V / p(r,r / ) 

s 

+ [n' x %¥(?')] x V'^r, r 7 )} dS' (2.73) 
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For 7 = 7*, the integrals in (2.72) and (2.73) are singular and by removing these singu- 
larity, they may be rewritten in terms of principal integrals as 

§■£(7) = T (7) + {-jfco [«' X rio'ff (7')] 5(7, 7*) 

+~ [ft' • V' x Tfoff (7 / )] V'5(7,T0 + [ft' x ’E(7 t )] x V^T')} dS f (2.74) 
l|jb7(7) = %7 i (7) + {jfco [»' X TCK)] 5(7, T') 

S' 

“ [n' • V x W)] V'5(7,n + [ft' x 1*7(7')] x V'5(7,7')| dS' (2.75) 

where we have also made use of (2.67) and (2.68). These must now be enforced on the 
boundary so that they can be coupled with the FEM equations. 

Initially, we will allow S' to be a general surface of revolution and will then speci aliz e it 
to the case of a right circular cylinder. In the next section, we derive the modal boundary 
integral equations by expressing the fields and the Green’s functions as a Fourier series 
in the cylindrical coordinate 4>. The resulting modal equations are then discretized and 
the resulting subsystem is augmented to the finite element system previously derived. 

2.2.1 Derivation of the Modal Boundary Integral Equation 

Consider the general surface of revolution indicated in fig. 2.1 whose tangential unit 
vectors are denoted by 4> and t. The angle v is that between the t and the z-axis and is 
negative when £ points toward the z-axis. Referring to the figure, we may represent the 
various unit vectors as 

ft = x cos v cos 4> + y cos v sin <f>— z sin v 
4> = — xsin^> + ycos<f> 


(2.76) 

(2.77) 
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£ = £ sin ucos^ + psinusin^ + f cost; (2.78) 

x = fsin uccw <£ + ncos vcos <£ - 'Asm <^> (2.79) 

y = i Bint; sin ^ + ncosvsin^+ 4>coe<f> (2.80) 

£ = {cost; — ft sin t; (2.81) 

Expressing the primed unit vectors in terms of the unprimed unit vectors results in 

i = i' [sin v' sin r cos(^ — <f/) + cos v cos t/] 

+n' [cos v' sin v cos (<f> — <j>') - cos v sin r'] + 4> [sin v sin(^ — ^)] (2.82) 

ft = i' [sin v' cos v cos(<f> — <f>') — sin t> cos v'] 

+«' [cos v 1 cos v cos(<f> — 4>') + sin v sin v'] + ft [cos v — ^)] (2.83) 

4> = —i' [sin v' sin(^> — ^)] — ft' [cos v' sin(<£ — <f/)} + 4> [cos(</> — 4>)} 

- -p' sin(<£ - <(>') + cos(^ - <f>) (2-84) 

Taking the <t> component of (2.74) and noting the identities, 

<j> • (n' x rjoH) = —TjoHjsmv' sin(4>— ft) - TjoH t co8(<f> — </>') (2.85) 

4> ■ Vg = -4> • Vg (2.86) 

n' • (V' x f h H) =j, [- UP'^ H 4) + £(*>**)] (2-87) 

j>- [(n'xEjx V'g] = 

[i'-Et sin(<£ - 4>') + n'Ej, cos (<f> - <f>') + 4/E4 cos v' sin(^ - ^')j • V'g (2.88) 
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we may rewrite (2.74) as 


5*W) = £j(F) 

/ / 0*0 [voE * Bin v' sin(</> - 4>) + tjqI T« cos(<£ - <j>)] g(T, f') 

J r J o 

+ ih [- ^' rioB * ) + ❖ww] * ■ v «( 7 - 7 ') 

+ [i'£t sin(^ - 4>) + h'E+ cos {<f> — 4>) + 4>E^ cos v'sin(^ - ^')J • V'j| p’dtfdT (2.89) 
Farther, by carrying out the derivatives of the Green ’6 functions, we have 

/ / { J*o [**> if# Bin tf' ainf^fr - <£') + t?o#i cos(^ - <j>')] g{7 t T f ) 

+ (Et sin(^ — 4> ') + E+W cos v' cos(<p — <f > r ) 

— p cos v' + (z — z') sin t>' cos(4> — ^')]) ~ I p'dj/dT (2.90) 

i£o J 

in which 

Ro - \[i? + p' 2 - 2pp* cos(<£ - ^) + (z - z*) 2 (2.91) 

To generate the corresponding integral equations for the modal components, the fields 
and Green’s function may be expanded as 


zoo= e 

ms— oo 

voS(t)= f; ^»(p,z)^ 

m=— oo 


n=— oo 


(2.92) 

(2.93) 

(2.94) 
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where 


?m(p,z) = j*_T:(p,u,z)e-i mu du 
Thn{p, z ) - 2% f Tf{p,u,z)t~ }mu du 


1 r* 

= 9n(p,p',z,z') - - f v co8(n«)dtt 
* J o 4*ri£ 

P^Obp'.*. 2 ') = | [Pn-l(P,P',*,*') + 0«+l(p,p', *,*')] 


1 f* e~^ R 

= — 4- cos « w- cos(nu)dti 

ir Jo 4?rfi 


0n 2) (P»P'>*> 2 ') = [Pn-l(p,p',2,*') ~ 0»+l(p,p', *,*')] 

= — + sin u =r- sin(nu)cfu 

it Jo 47ri2 

pi 0), (p>p', 2 , 2 ') = 0»(p,p\ 2 > 2 ') = ^2 / 0 ~ cos(»ti)du 

^^Cp^ 7 . 2 * 2 ') = 5 bn-lCP.P'. 2 . 27 ) + ^.+l(P»P , . 2 » 2 ')] 

1 r 1 dg . 

= 7 j f cos u-ar— v cos(nu )au 

^ o i2 <£iZ 

^'(P.P', 2 . 2 ') = [Pn-l(P,p'» 2 . 2 ') “ Pn+l(P» p , » 2 * 2 *)] 


-db/ 


r . i ^ . , x , 

U + sin «-=•—*? sm(nu)au 
q J o R dR 


R= yjp? + p n — 2ppf cos tt + (z — z') 2 


Substituting these into (2.90) yields 




+ E E *>*’*?’ + (‘ml9i 1, l 

n=-oom=-oo ^ 0 

[ — ^r(P 


(2.95) 

(2.96) 

(2.97) 

(2.98) 

(2.99) 
( 2 . 100 ) 

( 2 . 101 ) 

(2.102) 

(2.103) 
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-Cmt(z - ^)k^g^y + em+kltp' cosv'gW 
-p cos v'g n + (z - z') sin v'ffn ^)} p’d^'dT (2.104) 

and by multiplying each side by e~ ip * and integrating over (0,2*) to extract the mth 
modal equation results in 

2 e m4>(pi z ) — e \n<t,{Pi z ) 

+ 2 * j T {.X* sin v'gW - j ^(pX^' + G?X<) + rf'] 

-e mt {z - z')kogW + e m ^*o(p'cos v'pW’ 

-p cos v'g m + (z - z') sin v'p^')} fcop'rfT (2.105) 

after combining terms and where we have used 


/2ir I 

j d<j>' = | 


2* m = n 


0 otherwise 


(2.106) 


For the case in which T(= r o ) is the generating cross section of a right circular cylinder 
(indicated in fig. 2.2) the integral in (2.105) may be written as a sum of three integrals, 
one over each side of r„(= £?=i T^) as 


2 ®m^(p> z ) ^m^(p > z ) 

+ 2 * / 0 {j h m<t>g^} - j &r(p'hm4>)g$' + (jh mt ) [ p 2 ) +jmg^ 

-e mt (z - z 3 )k 0 ggy + e m4> k 0 (z - z 3 )p£>'} k 0 p'dp' 

+2T {" j ^bm^y + OX,) [«£> +jmgW] 

-e mt (z - z')kogW + e m *fc 0 (p2g$' ~ wl)} k 0 p 2 dz' 

+ 2 * / 0 {-jhm+gff - j ^(Am^Sm’ + C ?Xt) [s#* + jm*®'] 

-e mt (z - z')koggy - e m4> ko(z - z'JpW'} k 0 p'dp' (2.107) 
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Introducing the normalized coordinates 

R = kop R' = hop' 

Z = k 0 z Z ' = k 0 z' , (2.108) 

& = k 0 & 

(2.107) becomes 

±e m *(R,Z) = ei n t(R,Z) 

+ {jhm+gff - 3 ■&T(&hm*)g£ Y + (jhmt) [j^+rf'] 

+j(M(Z - Z 3 )gW + e m +{Z - Z 3 )g^'} R'dR' 

+ / 1 {-3 ■&r(R 2 f 'm*)g£ y + (jhmt) + jmg£ Y ] 

+iO«m.)(Z - Z')t£ y + - *ftj} *iZ' 

+ f* {-iAmwS 1 - i ,}>(.R'A m *)j£ ) ' + (jhmt) [si!* + jmsi?'] 

+ia«m.)(Z - Z')sS y - <W(Z - 7 , )sS ) '} KiR' (2.109) 

This equation and its dual are discretized in the next section. 

2.2.2 Discretization of the Modal Boundary Integral Equation 


Consider the fig. 2.2 where the rectangular boundary i6 divided into N a boundary 
elements and are equal in length along the r„ 2 . Along the boundary, the fields are 
expanded into pulse basis functions as 


N a 


£'(*'. z') = £ v j+i P(R i+i - R\ z )+! - z') 

:=i 

where U represents any one of the components e m ^, h m ^, e mi or h mt and 

1 iflZ j+ i-Z'\<^,\R j+i -R'\<± 


( 2 . 110 ) 


P(R jH -R\Z M -Z') = < 


( 2 . 111 ) 


0 otherwise 
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and 


Z i+i = i2i± V t2i R j+$ = (2.112) 

Substituting the pulse basis expansion into (2.109) and simplifying yields 

2 C *»*(-R> Z) — e m4>{ R i Z) 

+ E / {hm^} J+ i •£ _ jg^R'dR' 

j=N.i+N.3+l l J "i+i 

+{(R , h m 4 > )'} i+ ^ / - jg^ R'dR' + {jh mt }j -l* 3 [sm* + R’dR' 

J «>+i J 1 

HMif' /*' (z-z 3 ) iI wje'ije'} 

^ «>+» 3 *' Kj+i J 

JV«i+JV«a l z 

+ E ) {(*2>w)'} i+ i / - jgWRidZ’ 

j=N al + 1 l J Z i+i 

/*' [«£> + jmsfW] fl 2 dZ' 

J z J+1 1 J 

+{jc TOt }, - Z')gWR 2 dZ ' + {e m *} i+ £ - ^)iZ 2 ^| 

Wrt f r j* . . 

+ E | / - jg^R'dR’ 

i*= i l J 

+{(-R'h m ^)'} J+ i - jpW -R'd-R 7 + (ihmjj fl [jfi* + i^m^] 

p^ X o{Z - Z\ )g$' R'dR' - {e m *} j+ ^ f^*\z - Z x )gW R' dR^2.US) 
Proceeding to point-match at the boundary element midpoints, we have 


N. 


+ T. 

i=A r «i+A^«a+l 


Z^ ) — , Zi ) 

( {&!»*}*.$ jsOpR'dR' 


+{(*'/w)'} i+ i /*’ - mF-rw + { j/w}, /*■ f s m + jmp w| rn 

rti+l ./ itwo.1 1 J 
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+0«m.}»/ ,V j(Z; -Z 3 )gWR'dK' + {t„ i ) j+i P (Z,-Z 3 )sS}'R'dR'\ 

J H>+1 3 3 ^ ^41 3 ) 

{«*>**•<)'},+} -i£ y W 

+{jhmt}j i* kfi ) + R 2 dZ' 

J Zj + 1 1 1 

+{j*mth fz j+ J( z l - Z 'W R*dZ' + Mm ~ RmL)R*dz^ 

C * JL +1 

+ £ ] {&»»*},•+ 1 T - jg^R'dR' 

j=l l * J Rj 

+i(R^h m *) , }m /2 +l ■ j& ),R, dR’ + OW; /” >+1 [sS* + rf] fi'ifi' 

- Z x )gWR'dR! - {e m ,} i+| - Z{)g$' R'dR^llA) 

for the field points on contours r„i and and 

~e m ^(f£2j Z i+L ) = C m^(-^2) ^t+i- ) 

+ £ {{fcw} i+} /" < j&XiX 

j=H al +N*i+l K J A»+» 

+{(*X*)'} i+ j /*' - + O'w.Jj P l*S’ + wrg’l BOX 

+{;«-».}, P i(z,j - z^jftur + {«-*}** P (Z i+ i - Z3)jS>'J!'<i«'} 

*' «>+i * * J a J 

W«1 +JV02 /■ ,7- 

+ £ S {(-R2^m^)'}, + 1 f -jg$'R 2 dZ' 

:-Na 1+1 V ^ Z ->+» 

f ’ R*dZ' 

J Zj+1 1 j 

+{i'n.,h / Z '^j(z, +i - z')j£>'JMZ' + {e„*} )+J /;jfi JS w . 

+ £ i { W i+ * / - jg^R’dR! 

+{(f2 , /i m ^) , } J+ i J ^ - jgW&dR' + { jh mt }j [^ ) + -R'diZ' 
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+««-•}< /^‘j(Z, +i - Zi)l£ r B!dK - M i+i /"’*'(Z, +} - Z 1 )jW'fi'^115) 

for the held points on contour r a 2 « The above set of equations may be written more 
compactly as 

i [K -p*]c (^), - [/>■] - [i q*c + c*'b] {(w)„ - [g>] ovw,}. 

= { C m*} a | 2 . 116 ) 

where the matrix D arises from the derivative 

(2.117) 

and C is a matrix comprised of ’l’s along the diagonal and superdiagonal. Also, the 
subscript ^ represents evaluation at the boundary element midpoints. In a parallel 
fashion, the dual of corresponding (2.116) may be written 

i [W -P*]c {/w}„ - [f‘] {#W. + + q*'d ] + [g<] {,•««.}. 

The matrices in (2.116) or (2.118) are 3 x 3 in size, each element of which is a 
matrix corresponding a particular integration and observation (field point) contours. 
Each element of the submatrices is in Appendix F. For non-self-cell terms, the integrals 
are evaluated via open formula numerical integration schemes. The self-cell terms are 
given in detail in Appendix F. The integrals involving g m are computed via Romberg 
integration with a specified convergence criterion to ensure accurate evaluation for any 
mode. 
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Finally, augmenting finite element system with that fonned by (2.116) and (2.118), 


we derive the 6y6tem 


A U 

A il 

A ix 

A< Ia 

A U 

A ix 

A ‘xa 

A U 

A xx 

A da 

A h 

A dx 

P*C 

0 

0 

B aa 

B a i 

Bax 

Bu 

Bu 

Blx 

Bza 

Bu 

B zx 

Bda 

Bdi 

Bdx 

Q 

0 

0 


Kd 

Baa 

-Bai 

A )d 

-B Ia 

-Bn 

A U 

-Bza 

-Bu 

A h 

Bda 

-Bdi 

0 

-P X 

-Q 

B a d 

Ka 

<i 

Bid 

A 1a 

a ii 

B xd 

Ka 

A *u 

Bdd 

A da 

A di 


-Bax 

-Bad 

0 

-Bn 

-Bid 

0 

-Bzx 

-Bzd 

0 

-Bdx 

-Bdd 

0 

0 

0 

0 

A £z 

Ki 

Caa 

A ?x 

A Id 

0 

Kx 

A ?d 

0 

Kx 

A dd 

0 

0 

0 

0 


0 {<•*}.} 0 0 


c m^a 

n4>I 

jCmt 
hmfyi 

hm+I 

hmjd 
jhmt 

® ° {<*}.} ]( 2 U9 > 

( 2 . 120 ) 


0 Q* P*C 
= [ 0 0 0 

which is to be solved iteratively and where 

Q = Q+C + Q+'D 
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Chapter 3 

Scattered Field Computation 


In the far held the scattered fields are given by 


T(r) = E+(r)4>+rioHt(T)§ 
Vo~H‘(r) = rjoH^r)^ - £*(r)0 


We wish to compute the radar cross section given by [19] 

<r= H m 4x^Ell£=Uo,4x,>M 

For TM X polarization we have 


M 

E4>{r,0,4>) - 2j ^2 ej^(r, 6) sin(m^) 

msl 

Af 

ijo&}(r,0,<f>) = ho(r, 0) + 2 ^2 ft m*( r >^) c os(m^) 

m=l 

and for TE X polarization we have 

M 

E^{r,6,4>) = eo(r,0) + 2 ^ e^(r,tf)cos(m^) 

m=l 

M 

Tfc#J(r,0, 4>) = 2 j ^2 h m^( r y0)^ 

T7l=l 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 


28 



These combined with a unit amplitude incident field implies that (3.3) becomes 


’ M 2 M r 

= ton 2 53 «Ji*(r,0)sin(m^) + *$ + 2$: h^(r,$)cos(nuf>) (3.8) 

m= 1 m=l 

M 2 M 

cte m = ton 4xr 2 ej + 2 V e^(r, 0 ) cos(ro^) +2? h^(r, ^) sin(m^) 

~°° [ *-i 

We had previously discretized the Stratton-Chu integral equation for field points on 

the integration contour as given in (2.113). Eliminating the principle value factor for 

observation not on the integration contour, the corresponding scattered field equation 

may be written 

e' m4> (R,Z) = j {jKn+gW - j + (jh mt ) [$£> + jmg£ y ] 

+j(je m t)(Z - Z 3 )gW + e m +(Z - Z a )g£Y} R'dR ' 

+ / {-J ■&t(Ri h m+)9 { m Y + (jh m i ) WP + 

+i(fc*X* - z'lsiS’' + - *»!.)} W 

+ -i + (jhmt) [sS* + jrojji?'] 

+j(;«m,)(Z - z'JsS 1, - «m*(2 - Z')si! r } •R'd.R' (3.10) 

We wish to evaluate this expression for large kor = VR 2 + Z 2 . For large r 

/ 

Jr 2 + £' 2 - 2fLR' cos u + (Z - Z') 2 ~ k 0 r - =±- - cos u (3.11) 

v kor kor 

Thus, we may write (2.97) 

g m (R,R\Z,Z') as — — e J &‘ [* c , ^' co * u cos (mu)du (3.12) 

jt 2*or Jo 

Noting that the integral is related to the Bessel function of the first kind, we may write 
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(3.12) as 


g m (R,R',Z,Z') ~ -^j—eiZ'^'fmiR'J) 

where we have used 

= ~ r e* 0aMX cos{mx)dx 
7T J o 

= MU*) 

and the fact that 


JZ = kor sin 9 Z = Jfc 0 r cos 0 
Likewise (2.98) - (2.102) become 

a W(R,K,Z,Z') at £0e>' z ' 

sL I| (A,-R',z,z') a £0^—'/™^,*) 
sUR.B'iZ.Z') a -~^0 e ' 2 ’“V m («',9) 
j£>'(*-R',Z,Z') s 

where 


fcm(R?,$) = /""^(iZ'sintf) 

where the prime on J indicates differentiation with respect to the argument, 
these expressions into (3.10) results in the expression 

g-ifcor 

c m$(-^>^) = 2k^r 


(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 

(3.23) 

Substituting 

(3.24) 
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where 


ih m +(jf m ) + (jhmt)fcm + (j^mt) CO B0f m 

Jo 

-Cm* COB 0(jfem)} R' dR' 

+ f* * Z ' C - 9 {(jk m t)fcm + (JCmi) COS 0fim + ^ sin*(j/ m )} R 2 dZ’ 

JZa 

+e jZt ccl {-hn+VM + {jhmt)fcm + (J^t)cOB0f m 

Jo 

+e m <t> coe 0(jfcm)} R'dR' 

Using a midpoint integration to compute the integrals, (3.25) becomes 

fe(m,0) = e? 2 *™ 9 £ {{h m +} j+ i\jf. m (R i+ i,0)} + iMjfcm 

j=N.i+AT. a+1 

COS Oftm 1 COS 0(J fcm ) ^ Rj+L Aj 

N+i+Nca 

+ 53 * ’'t** {{j h mt}jfcm(R2,0)+{je m t}jCOe9f m 

jsTV.j+1 

+{<w}i+i rintf(i/ m )} iZ a Aj 

w«, 

1 53 \~{^m^}y+lL7"/»»i»(-Rj+l*^)] + O^rafJj/em + 0" e ml}i COB0f tm 

i« l 

Letting A(m,0) be the dual of (3.25) we may write (3.8) as 

2 E r/,(m,«)sin(m«l 


Aq 4jt 


m=l 

M 


< T TE.( 9 y 4 >) 1 _ 

Ag 4 X 


AM + 2 £ i m A(m, 0) cos(m^) 

m=l 

M 

fe(0,9) + 2 53 J m /e(m, 0) cos(m^) 
m=l 

M 

2 53 i*" A( m » 0) sin(m^) 


m=l 


(3.25) 


(3.26) 


(3.27) 


(3.28) 
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where Aj is the length of the jth boundary element. 
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Chapter 4 

Results 


The scattering patterns for a test body are shown in figs. 4.1 - 4.4. The structure 
is a conducting right circular cylinder of length 1A and radius ^j. Fig. 4.1 shows both 
the TE and TM cases for broadside incidence for mode 0 and as seen these are in good 
agreement with corresponding data based on the MOM code CICERO [7], except for 
some deviation of the TE curve in the region between 0 and 30 degrees. The results for 
mode 1 are shown in fig. 4.2 and again, similar observations are applicable in this case 
as well. Fig. 4.3 shows the sum of modes 0 and 1, where we now observe a disagreement 
of the TM curves indicating that the phase associated with modes 0 and 1 must differ 
with respect to the data obtained from the CICERO code. 

Fig. 4.4 shows the bistatic scattering pattern for the Bame geometry with axial 
incidence. Only mode 1 yields a non-zero solution in this case and the depicted results 
again show some deviation from the reference data in the forward and backscattering 
regions. Presently, we are investigating the cause of these disagreements and in addition, 
we are researching new approaches to improve the storage and computational efficiency 
of our code. 
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ctA [dB] 









Figure 4.1: Mode 0 TM and TE bistatic scattering pattern from a perfectly conducting 
circular cylinder of length 1A and radius O.lA for broadside incidence. 
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oA [dB] 


p=.lX, 1=1. OX pc Cylinder (m=l) 



6, [deg] (^=0,8-90) 


FE/BE (TE) 

FE/BE (TM) 

© CICERO (TE) 

a CICERO (TM) 


Figure 4.2: Mode 1 TAf and TE bistatic scattering pattern from a perfectly conducting 
circular cylinder of length 1A and radius 0.1A for broadside incidence. 
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lap] r xA> 


p=.lX, 1=1. OX pc Cylinder (m=0+l) 



Figure 4.3: Modes 0+1 TM and TE bistatic scattering pattern from a perfectly con' 
ducting circular cylinder of length 1A and radius 0.1A for broadside incidence. 
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aA [dB] 


p=.lX, 1=1. OX pc Cylinder 



0, [deg] (^=0,6=0) 


FE/BE (TE) 

— FE/BE (TM) 

© CICERO (TE) 

b CICERO (TM) 


Figure 4.4: TM and TE bistatic scattering pattern from a perfectly conducting circular 
cylinder of length 1A and radius 0.1A for axial incidence. 
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Appendix A 

Derivation of Modal Incident 
Field 


Consider a field incident at a point T = (r,<j>,z) at an angle (O', ft) (as indicated in 
fig. 2.1 ) of the form 

(A.l) 

?(*•', ft; p, ft z) = -Ve’*** (A.2) 

where the ft direction is perpendicular to the plane of incidence and O' direction is in 
the plane of incidence. Using 

f = x sin# cos ^ + ysin0sin^ + icosl (A.3) 

z! = x sin O' cos ft + y sin O' sin ft + z cos O' ( A.4) 

the argument of the exponential becomes 

F 0 .f = k 0 r(-z’-r) (A. 5) 

as — fco [psinfl' cos(^ — ft) + zcosfl'j (A.6) 
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in cylindrical coordinate system. Using these and the fact that 


= — £sin0* + ycosfl* 

(A.7) 

= f cos 0* cos <f>' + y cos 0* sin <f>* - z sin O' 

(A.8) 

A A 

f = f sin 0 cos 4> + 0 co6 0 cos $ sin 

(A.9) 

y = f sin 0 sin ^ + 0 cos 0 sin <f> + ^cos ^ 

(A.10) 

i — r cos 0 — 0 sin 0 

(A.11) 


(A.l) and (A.2) become 

4> - <t>\ z) = [psm(4> - 4?) + 4>cos{<t> - *•')] co.fi>] (A 12) 

£{0'\p, <t> — 4>\ z) = — Jp cos O' cos(^ — <f > ') — ^cos O' sin(^ — <f> % ) — Z sin 0’j 

e >^[/>dnfi>co.(^— ^*)+x co.fi>] 


The previously derived fields may be expanded into a Fourier series in the parameter 


4> — <f> 1 by first writing (A.l) and (A.2) as 

?(*w -*’,*)= e (a. 14) 

me-oo 

W-,p,t-4,\z)= f. (A.15) 

m=— oo 

and then making the definitions 

/(**; p, <t> - 4?) = “ <*-*) (A. 16) 

fc(0';p,4>- V) = cos(<£ - ^’)/(^;p,^- <P) (A. 17) 

/.(«’'; P,4> ~ V) = «n(* - Wf^pA- **') (A. 18) 
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Expanding each of these into a Fourier series in (<f> — 4> x ) and using the fact that 



M) = /(-*) ^ /m(u) = /- m («) 

(A.19) 


f{4>) = <=* fm(u) = -Mu) 

(A.20) 

we have 

<t> { ) = M»\p) + 2 2 fm(0 i ,p)coB[m(<l> - 4>*)] 

m=l 

(A.21) 


M&’iPyt- <f>') = fd0(6\p) + 2 53 fcm(9\p)cOB[m(<i> ~ ^‘)] 

m=l 

(A.22) 


UV\P,4>- f) = 2 j £ /, ro (0)sin[m(^ - **')] 

m=l 

(A.23) 

where 

f m (6\p) = - /V*'“*~“cos(mu)du 
* Jo 

(A.24) 


fcm(B\p) = - r cos co,, ‘ cos(mu)du 

IT Jo 

(A.25) 


f m (6\p) = [* sin uc j «>•« gin(mu)du 

5T Jo 

(A.26) 

Noting the identities 



f 71 Jm(fl) = — f ccm x cos(mr)di 

JT JO 

(A.27) 


j m_1 = — / cob cos(mz)dx 

Jo 

(A.28) 


C sin xe^ 00 ** sin(mi)di 

(A.29) 

where the last two are derived from the first by differentiation with respect to /3 and 

integration by parts respectively, we may write (A.24)-(A.26) as 



fm{0\p) ~ f 1 JmiksPSluO') 

(A.30) 
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With these, we may proceed to rewrite (A.14) as 

m=— oo 

£ [p/cm(0\p) coed 1 - famitf'p) cos $ — i/m(0*,p)sin 

ms— oo 

or, using (A. 19) and (A.20), we have 

[/52j /»m(0V)sin[m(^- 


+4>fdo(6 i ,p) + <^2 23 /cm(^‘,p) cos[m(^ - 


msl 
oo 


msl 


•|pCO 6 0 ‘ 


— isin0* 


fto(0\p) + 2 2Z /cm(^',p)cos[m(^-d>‘')] 


~4>co& O' 


msl 

oo 


2 j J3 Am(0 f »P) sin[m(^ - 4>*)] 


m=l 

oo 


fo(6\p) + 2 23 /m(0‘,p)cos[m(<£ 


m*=l 


-*•))]} 


In this work, we will use the <f> components of each of these equations. 


(A.31) 

(A.32) 

(A.33) 

(A.34) 

(A.35) 

(A.36) 
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Appendix B 


Maxwell’s Equations for 
Axisymmetric Media 


The usual Maxwell’s equations in a source free region are given by 


V x £(r) = -jupH (B.l) 

V x ~H(r) = jueE (B.2) 

V.E(7) = 0 (B.3) 

V-'5(7) = 0 (B.4) 


In cylindrical coordinates the electric and magnetic fields may be expanded into a Fourier 
series in <f> as 


E (b.s) 

ms-oo 

vE(r)= 52 Hm(p,zy m4> (B.6) 

ms-oo 

Substituting these into Maxwell’s equations, we obtain 

Vrt ,- 3 —? m x i = 

P V 

VxL + —h m X <£ = ju)€TjZ n 

P 
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(B.7) 

(B.8) 



Appendix C 

Derivation of Boundary 
Conditions 


In this appendix, the axial and perfectly conducting boundary conditions are derived. 

C.l Derivation of Axial Boundary Conditions 

Substituting the Fourier series representation of the electric field into the divergence 
condition we obtain in the normalized cylindrical coordinate system 

V-(? m ei' n *) = k 0 e’ m +{j l e r e mp + &(e r 

^mp ) + j%e r e m 4 > + &( 

^r^mi)] — 0 (c.i) 

Thus, 

€ r e mp + R [ &(e r Cmp) + ig(*rt m z)] = -jTn€ r e m ^ (C.2) 

as Morgan had previously derived. Talcing the limit of this expression as R — ► 0 + , we 
obtain 

Cmp + jme m * = 0 (C.3) 

Expanding the derivative w.r.t. R in 

i U me mp = Jtbh mz (C- 4) 
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and taking the limit as R — * 0 + , we obtain 


®mp — 0 (C.5) 

Combining (C.3) and (C.5) and solving for e m ^ we have 

(m 2 - = 0 (C.6) 

In a similar manner we obtain the dual expression 

(to 2 - l)/WI*-o+ = 0 (C.7) 

For m ^ 1, the following axial condition holds 

|ft=0+ = ^m^l/i=o+ = 0 (fft 7^ 1) (C.8) 

To derive the condition for m = 1, lets first consider 

e mz = jfm & (Be m ^) “ B/i r (B/l m ^)] (C.9) 


As B — ► 0 + , e mx — 1 ► 0 for m 0. Differentiating (C.9) with respect to Z we have 
^T c mx — jfm — B [ + B ^^m^) 

+/*r ( + B } (C.10) 

Clearly, as R — ♦ 0 + -fetmi = 0 for m / 0. Differentiating (C.2) with respect to R after 
dividing by e r , we obtain 

"ik e mp +£ [ ^r^mp) + |/fe=0+ = ^^m^l/lsO+ (C.ll) 

Accounting for the behavior of e mx and (C.ll) becomes 

2 “1" e mp &(ln c r ) + ;m &e m * = 0 (C.12) 
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To find another equation in terms of ^e mp and we multiply (C.4) by R 


and differentiate it with respect to R to obtain 

R A0* r ^”**) + Hrhmz = j [2 + R- ^T e m^ j m A c "»p] (C.13) 

Letting R-+ 0 + we obtain 

2 j AW + m -§£tm4> = 0 (C.14) 

Substituting (C.14) and (C.5) into (C.12) we obtain 

(4-™ 2 ) Aw + W A(her) = 0 R = 0, m # 0 (C.15) 

In an analogous fashion, the dual of (C.15) is given by 

( 4 ~ m2 ) A*W + ^(In/ir) = 0 R = 0, m f 0 (C.16) 

For c r and fi r constant in R at the axis of symmetry and for m = 1, (C.15) and (C.16) 
reduce to 

AW = 0 (C.17) 

A*W = 0 (C.18) 

C.2 Derivation of PEC Boundary Conditions 

On a perfect conductor the condition 

nxE= 0 (C.19) 

Substituting the Fourier series expansion for the field into this boundary condition yields 
the following condition on each mode 

nx? m =0 (C.20) 


48 



(C.21) 


The second Maxwell’s equation for the mth mode is given by (see an appendix) 

V x Tl m t hm x i — 

P 

Crossing this equation with n and noting that n • hm = 0 on the conducting surface, we 
obtain 

n x (V x Tim) — 0 (C.22) 

Carrying out the curl in cylindrical coordates yields 

ft * {p ^ ~ 

+* [ ^(P^m*) — p} = ® (C.23) 

Noting the identities 


n X p = 4>{h • z) 

(C.24) 

n x 4> = * 

(C.25) 

n x z = — $(n • p) 

(C.26) 


we And that the middle term of (C.23) implies 

&tlmp — ‘fahm* (C.27) 

and the first and third terms may be written 

^(6 ' *) [p "£fhmx ~ ~ 4>{ft • p) | 'fciphmt) ~ * = 0 (C.28) 

Rearranging terms, we have 

ft • [V t (ph m +) - 1 £h mt ] = 0 (C.29) 
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or 


= o 


(C.30) 


and we have used the following 


= + (C.31) 

£ = ft*V« (0.32) 

V’/i = koph m <j> (C.33) 
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Appendix D 

Evaluation of Finite Element 
Contour Integral 

D.l Contour Integral Evaluation along Conducting Sur- 
faces 

It is shown in the appendix that along perfectly conducting surfaces the conditions 

V»« = 0 (D.l) 

= 0 (D.2) 

must hold. During the assembly of the finite element equations (i.e., when the summation 
over all elements is performed), those rows and columns of the finite element matrix 
corresponding to nodes on the conducting boundary are eliminated. As a result, the 
corresponding contour integral vanishes along a conducting boundary. 

Imposing the condition (D.2) results in the elimination of the associated contour 
integral since on the conducting surface 

n • (£ x Vt&) = t • V t ^ e = 0 (D.3) 

(check this stuff) 
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D.2 Contour Integral Evaluation along the Axis of Sym- 
metry 

In the appendix, the axial boundary conditions are derived 


«w * 0 

(D.4) 

o 

II 

J 

(D.5) 

= 0 

(D.6) 

— 0 

(D.7) 


Conditions (D.4) and (D.5) results in the elimination of the rows and columns of the 
assembled finite element matrix associated with nodes on the axis. 

Alternatively, since R -► 0 all terms in the contour integral axe zero by virtue of the 
chosen weighting function. 

(may explore the possiblity of a different weighting function which does not guarantee 
this) 

D.3 Contour Integral Inter-element Connection Cancel- 
lation 

Since the argument of the contour integrals are tangential fields at the element bound- 
ary, they will be continuous between adjacent elements. As a result, the contour inte- 
grations along the element intersection cancel. 
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Appendix E 

Evaluation of the Finite Element 
Matrix Elements 


In the evaluation of a- and b'j, integrals of the form 


and 


Pab ~ 



R a Z b dRdZ 


(E.1) 



R a Z h 

R*k 2 - m 2 


dRdZ 


(E.2) 


Clearly, Q ab exhibits singularities for real k. To evaluate this integral, consider an integral 
of the form 


1= JJ g(R)Z b dRdZ (E.3) 

S' 

To evaluate the integral, first transform it to an integration along the element boundary 
via 


Using Stokes’ theorem 


zH-i 

v(R,Z) = -g(R,Z)f^jp 


JJ ( V x v) • dS = / T) • dl 

n » C* 
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(E.4) 


(E-5) 


f 



and 


Vx17= -j>g(R)Z b 
dS = 4>dRdZ 


(E.6) 

(E.7) 


Inserting these into (E.5) yields 

J j g(R)Z b dRdZ = j> g(R)Z b+1 dR (E.8) 

S* ^ 

Via (E.8), (E.l) and (E.2) become respectively 

Pat = t-J-t / R a Z b+1 dR (E-9) 

o + 1 Jr* 

and 

Qab = — l f a Z b -?—dR (E.10) 

Voi 6+1 Jr* R 2 k 2 -m 2 K ’ 

where the contour integration is taken in a counterclockwise fashion. For linear triangular 

elements, T e is represented by a summation of three contours, one for each side of the 

triangle. The variable Z may be thus expressed as 

Z(R) = u t R + vi (E.ll) 


where 

_ Z/+i - Z\ 

~ R, +l - Ri 

v\ = Z\ — tuRi 

Then Z^ 1 may be expressed as 

Z* 1 = (“+** >) WI = + 



(E.12) 

(E.13) 


(E.14) 
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where 


h- ” ! 

\m/ m!(n — m)! 


(E.15) 


Thus, by writing the integral in (E.10) as a sum of an integral along each side of the 
triangular element, it be rewritten as 

3 3+1 r: 1 ) ft )>r 

. psO 

Clearly, integrals of the form 




R 2 k 2 — m 2 


dJt 


(E.16) 


/(n ’ m) -/ R , 


# 2 * 2 - m 2 

for n = 0, 1, ..., a + b + 1 must be solved. For n = 0 

l^+i 

I 

J(0, m) = 


dR 


(E.17) 


Hi 


2m»t P®( m “ **) ~ ln ( m + -R«)] 

For n s 1 it is easily shown that 

7(1, m) = [ln(»n - Rk) + ln(m + Jte)] 


fk+i 


m = 0 


m > 0 


(E.18) 


Hi+i 


m > 0 


(E.19) 


Using the definition of the principle branch of the natural logarithm in the equations 
above guarantees that the singularity is properly handled. For values of n > 1, the 
recursive formula [17] 

,2 


m 


7(n, m) = 7(n, 0) + — =-7(n - 2, m) 


(E.20) 


is used. Thus, (E.16) may be written in terms of 7 as 

Ttgo; 1 ) (*)>'<»+>>■-»> 


(E.21) 
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In a similar fashion, 


The shape function is written in expanded form as 


Hi+i 

* 


where 


(E.22) 


(E.23) 


We had derived in section 2.1 


n* = 

(E.24) 

o' = ZSK-ZiX 

(E.25) 

•*<» 

II 

1 

& 

(E.26) 

7 ? = 22-3 

(E.27) 

RV t (RN?) ■ V t (RNf) + c r iV?JV/] dS* 

(E.28) 


Noting that 

ff'N! = jjjfp [ofoj + w°’j + + Xl! aj + 1j°!) 

+R2(Ph; + d’pt) + z 3 PtPj + flS.'-i, - ] 

Substituting these into (E.28) and reducing we obtain the desired result 

«j h = + Wu - 2(7*0* + 7 ;O*)Q20 

-2(^7* + Pj7i)Q 21 - W'Ql2 - (47*7' + /?*)<?30 


(E.29) 


(E.30) 
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(E.31) 


+«?«; *0 + (#<*' + FiCt\)Pi 1 + (7f«I + 7;«f)P20 + (tf 7, ? + m !)Pn 

+PtPjP 13 + 7?7i-P3o] (2 ^ e)2 

In a similar maimer, we may write 

b v * II M x • W Mi] 


as 




and likewise as 


*& = + HP*i* - PhDQ»\ 


(E.32) 


(E.33) 


(E.34) 
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Appendix F 


Boundary Integral Matrix 
Elements 


In this appendix, the elements for the discrete boundary integral system axe pre- 
sented. 


F.l Elements of P* 


ft]. 

= 0 

(F.l) 

hi 

= f’ [ft -Rn. i s'M H ,R I ,Z 1 ,Z')}R^Z' 

(F.2) 


= (Zj - Z 3 ) 4** g£Y(R i+ L,R',Z l ,Z3)R'dR' 

J fy+i * 

(F.3) 

Hi, 

rRj+i 

= (z^j-z, )/ R 

(F.4) 

H}* 

= f* [jStfg*' - Z i+ j,Z')] JWT 

(F.5) 

Hi, 

= (Z j+i - Z,) /*’ «g ) '(K 2 ,fl’,Z j+ i,Z3)J!'<U!' 

(F.6) 

Kl, 

= (Z 3 - ZO /^ +I 

(F.7) 
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(F.8) 


[/?,].. = f‘ - «. +} s;(«, + J,Fl,Z3,Z')] R,iZ' 

= 0 <*•») 

F.2 Elements of F* 


[«*]„ - 

0 

(F.10) 

N, - 

f Zj m - Z')gW(R i+i ,R 2 ,Zi,Z')R 2 dZ' 

J Zj+ 1 

(P.11) 

Cl 

ll 

(Zi ~ Z 3 ) (** jgW{R i+ i,R\Z l ,Z 3 )R'dR' 

(F.12) 

w, = 

(Z i+ 1 - Zi)J^ 39$' (Rit R'l Z i+ 1 , Z\)R!dR! 

(F.13) 

[ p »],i = 

C i( z i+x - Z')g£Y(R 2 ,R 2 ,Z i+k ,Z')R 2 dZ' 
J Zj + 1 2 3 

(F.14) 

[ P 23 L; = 

(Z i+ 1 - 23 ) /* tf,Z j+ ,,Z, > )*'<UJ' 

2 JR j + 1 3 

(F.15) 

M, = 

(Z 3 - Zi) /^ +1 jgW(R i ^R l ,Z 3 ,Z 1 )R , dR f 

(F.16) 

[ P32 ]y = 

i(Z 3 - Z')gW(R i+ ± , R 2 , Z 3 , Z')* 2 dZ' 

(F.17) 

l- o 1 
C: 

II 

0 

(F.18) 

Elements of Q * 


[«»]« 

= /^ +1 - jg% ) (R t+h ,R',Z u Z 1 )R’dR > 

J Rj 3 

(F.19) 

[«-],, 

= 0 

(F.20) 

[««],, 

= f**'jg% ) (R i+ L , Z u Z 3 )R’dR' 

(F.21) 
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[(&].. = (F.22) 

[<&]„ = 0 (*■»> 

[<&].. = X.Z^ZiiKiK (F.24) 

[gl,].. = -«£>(«,>}, Jf', a, Zl)^ (F-25) 

[<?»],.,. = 0 (F.26) 

[<&].. - l^'jsSKR^.R'.Zs.Z^dR’ (F.27) 

^ ^ U J Rj 

F.4 Elements of Q *' 


[«?;], 

= 

/*' +l -39% y (R i+ ^R',ZuZi)R'dR' 
J Rj 

(F.28) 

[««]« 

= 

f‘ -i9 l £ f i.R, + i,R^z l ,z')Z'iZ' 
J Z«i 1 

(F.29) 

K] (i 

= 

f‘" - jsg> , (*, + i,a',Zi,z 3 )*'<«' 

J Rj 

(F.30) 

[«S]« 

= 

■[**' -jg^\R 2 ,R , ,Z iH ,Z 1 )R'dR l 

J Rj 

(F.31) 

[ Q ^]a 

s 

f’ -jgWiR^R^Z^^Z'V'dZ' 

(F.32) 

[<*]„ 

= 

Z^ +l -j 5 L 2) '(^,-R',^ + i,^3)iiW 

•/ fl> a 

(F.33) 


35 


(F.34) 

[®»]« 

= 

j Zj+1 2 

(F.35) 

[«£]„ 

5= 

/* ,+1 -i^'(fl, + i,3R , ,Z 3> Z 3 )iiW 

(F.36) 
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F.5 Elements of Q * 


[<?u] 0 . = /^ +1 [a™ +i m & ) \R i +L,R\ZuZ l )] R'dR' (F.37) 

[$«]- = (F.38) 

[Q 13 ] 0 - = /^ +1 + j m 9m Y ( R i+$ * R, y Z i > Z 3 j\ R'dR* (F.39) 

[Q 11 ] f . = ^ 1 [9^ ) +jmgW(R 2 ,R r ,Z i+i ,Z 1 )]lt l dff (F.40) 

[<&].. = {*’ {sW + jTng£ Y (R2,R2,Z i+k ,Z')] Z'dZ' (F.41) 

[<?»] u = /^ +1 [*£ } + jmgW(R 2 , R', Z i+ i , Z 3 ) ] R'dR! (F.42) 

[di] 0 - - j ^ [pL 1} + j™9 ( m Y (R i+ L , R', z 3 , Zt )] R'dR' (F.43) 

[Qla].. = [aW+jrngWiR^R^Z^Z')} Z'dZ' (F.44) 

v J Zj+i 1 2 J 

[^ 33 ] 0 - = /^ +1 l^m+Mg^ y (R i+ x,R',Z 3 ,Z 3 )J R’dR' (F.45) 

F.6 Self-Cell Evaluation 


The integrals in the matrix elements [P*],„ [<?£],-, and [Qy,-,- contain integrable 
singularities. They could be integrated numerically without modification as long as the 
singularity point is avoided, but costs excessive computation time. To avoid the resulting 
excessive computation time and innacuracies, the integrals are evaluated as in [5]. 

For self-cell integrals involving g m , Glisson gives 

th 

/, , R', Z i+ i , Z')R'dl’ 

*/ ij * * 
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(F.46) 


where 


= - /' f 

2jt J /, J o 


e~iR° . D , 

— = — cosfmtijiT 

Ro Ro 


dudl' 

+I(R i+ ^,l i+ i,luh) 


r/n 1 1 . , j>+l [ h [ 1 jr / , i^j_] 

i(R,i,h,h)- w t h [ n7 K y n 2 ) 2ii i+ | 

+r- l(*a - M - (la - Wfa -*)-(*- MM* - *i)] 


.(P.47) 


and where jiT is the complete elliptical integral of the fast kind, I may be either Z ot R 
and 


78, = ,/(*.+} - W + (Z i+ J " Z ')’ (F-48) 

78a = V’C-R.+i + *') ! + (««.} - Z’) 1 ( F - 4 9) 

Also, 

Ifo = + -R' 2 - 2J?- + iiZcoBu + (Z i+ ^ - Z'f (F.50) 

(F.51) 


The first and second integrals of (F.46) may be computed using an open interval numer- 
ical scheme that also avoids the midpoint. 

The integral expression for the self-cell of P*? may be rewritten as 


(l + jflo)e cosim^du 

2Rq 


[ p -],= /z‘ +1 [--Jr su % 

i r ( i+jr ' 
x 7 o R$ 

= I Z< - sin 2 (^)cos(mtf)rfu 

7 z, + i x 7 o R% Ro 2 


1 r(i+j|el£^l cos(m „ )d „ 

2Rq 


R%dZ' 

R\dZ' 


(F.52) 
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where we have used the identity 


1 — cos u = 2 sin 2 (^ 


The solution to (F.52) is 


i / 

•* J Zi+J * J i 


(1 + jRo) e~*^° , N . 2 , 

— g ~ " ' fe cos(mu) sin 2 ($ 


Ro 

.2 


tr 


4 + (Z, +i - Z')i] 


3/2 


duB*dZ' 


+&r(R 2 ,z i+ i,z i+l ,Zi) 


where 


I'W,luh) = /’ f 

J h J o 


:dudV 


hJo [iZ2„2 + (/ _ l')2 ] 3 /2 

= {(/ - /x)ln [ifcr + yj(l - 1 1 )2 + ^2,r2j 

+(/ 3 - /)ln [ifcr + v /(/ 2 -/) 2 + £ 2 5 r 2 ] 

-(/ - /i)ln(/ - h) - (/ 2 - /)ln(/ 2 - 0} 


In the same manner we have 




(1 +jRo)e->*> . , 

-b — sin(mu) sin(u)J2 


R% 2Ro 

™ 2 R i+ i 


duRfdl' 


2 [^« 2 + (Je,. +i - J R , )2] 


(F.53) 


(F.54) 


(F.55) 


(F.56) 
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where 


h = R< 

h = iZj+i 

li+$ = R i+\ 

V = R' 

for a = 1 

h = Z«+i 

h = Zi 

*•+§ = Z i+$ 

V = Z' 

for a s 2 

h — R \+ 1 

h - Ri 

f t+ ^ = -R,- + £ 

r = r' 

for a = 3 


Finally, we treat each term in [<?£*],,• seperately and obtain 


[<?»«],.. = [rf? + J' m 5rn ) ' (R i+ $ > #1 Z i+± » Z ')} R ' dR ' 

= J- I h f \l^1 cob tt cos(mu)R' - dudV 

2* Jh Jo Ro Ho 


r (l+JR0)e~ JR0 . , \ . jyf 

- - — sr — - — =: — sin(mu) sin uR 

0 R% 2Rq 


mu 2 R- 


duRTdl' 




(F.57) 


(F.58) 


where (F.57) is used to determine the expression for each value of a. 

The self cells involved in the other matrices contain non-singular integrands and may 
thus be integrated numerically without modification. 
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